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Abstract 

Among the direct search experiments for weakly interacting massive parti¬ 
cle (WIMP) dark matter, the DAMA experiment observed an annual mod¬ 
ulation signal interpreted as WIMP interactions with a significance of 9.2cr. 
Recently, Jonathan Davis claimed that the DAMA modulation may be in¬ 
terpreted on the basis of the neutron scattering events induced by the muons 
and neutrinos together. We tried to simulate the neutron backgrounds at 
the Gran Sasso and Yangyang laboratory with and without the polyethy¬ 
lene shielding to quantify the effects of the ambient neutrons on the direct 
detection experiments based on the crystals. 


1. Introduction 


Numerous astronomical observations have led to the conclusion that the 
majority of the matter in our universe is invisible, exotic, and nonrelativistic 
dark matter 0,3- However, it is still unknown what the dark matter is. 
Weakly interacting massive particles (WIMPs) are one of the most attractive 
dark matter particle candidates [3;, 0]. The lightest supersymmetric particle 
(LSP) is a possible WIMP candidate theoretically hypothesized beyond the 
standard model of particle physics. There have been a number of experiments 
that directly search for WIMPs in our galaxy by looking for WIMP-nucleus 
scattering through nuclear recoil 0 :6]. 

To date, several experiments (e.g., DAMA 0,3) CoGeNT [12 Jjjj], CRESST 
0 , and CDMS 0) have yielded signals that could be interpreted as being 
possibly due to WIMP interactions. Except for the DAMA experiment, the 
indications of WIMP-like signals are not significant at the <3a confidence 
level, and the same experimental groups or independent groups are planning 
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to confirm these signals. However, the result from the DAMA experiment 
has attracted attention because the observation of an annual modulation of 
WIMP-like signals with a significance of 9.2a has been reported. This find¬ 
ing has spurred a continuing debate concerning the observation of WIMPs 
over the past 15 years. The WIMP-nucleon cross sections inferred from the 
DAMA modulation are in conflict with limits from other experiments that di¬ 
rectly measure the nuclear recoil signals, such as XENONIOO 17], LUX 18], 
and SuperCDMS [jj| . However, it is possible to explain all these experi¬ 
mental results without conflict because of nontrivial systematic differences 


the WIMP distribution 22 


in detector responses 20, 21] and the commonly used astronomical model for 


Recently Jonathan Davis of Durham University has proposed a new model 
for the DAMA annual modulation, which is a sum of two annually modu¬ 
lating components with different phases. More specifically, the events are 
composed of neutrons, which are liberated in the material surrounding the 
detector by a combination of 8B solar neutrinos and atmospheric muons [?]. 
The DAMA group disagreed with Davis’s claim because the induced modu¬ 
lation amplitudes from neutrons induced by muons and by neutrinos are less 
than 9 xlO” 6 cpd/kg/keV (2 xlO” 5 cpd/kg/keV for neutrons produced in 
the lead shield) and less than 2 xlO” 6 cpd/kg/keV, respectively from the 
simulation, and they are less than 0.1% of the modulation amplitude mea¬ 
sured by DAMA/LIBRA 0- Further, Barbeau et al criticized Davis’s claim 
because a seven order of magnitude discrepancy in the neutron contribution 
was required 


2. Method of simulation 

To understand the neutron backgrounds for Nal dark matter experiments, 
we have performed simulations with the GEANT4 Toolkit 37]. The hadronic 


models of muon nucleus interaction and muon decay at rest were used in 
the GEANT4 simulation, version of 9.6.p02 (G4MuonVDNuclearModel and 
G4MuonMinusCaptureAtRest). 

We have simulated the cosmic muons and the neutrons induced by sim¬ 
ulating muons passing through the rock directly. We did not simulate the 
scintillation and light collection processes inside the crystals because it is 
not required for the conclusion derived in this study. Further, we did not 
simulate the environment neutrons underground because presumably these 
neutrons will not have annual modulation signals. 
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3. Analysis 


3.1. Muon energy spectrum 

Neutrons are produced primarily by the muons in the rock and materials 
in the detection system. These neutrons generate secondary neutrons via 
hadronic interactions in the materials. We simulated all the primary and 
secondary processes using GEANT4 simulation program starting with muons 
passing through the rock. 


We used the muon energy spectrum provided by D. M. Mei and A. Hime 34 
which is given by 


dN 

dE,, 


= Ae- bh ^- 1 ^{E il + e^l-e~ bh ))- 


■7m 


( 1 ) 


where h is the rock slant depth in km.w.e and we used a set of parameters 
provided by Groom et al. for e M , b, and 7 ^. 



Figure 1: Muon energy spectrum at Gran Sasso (blue) and at Yangyang (red). 


For two different depths at Gran Sasso and at Yangyang laboratory (Y2L) 
we assumed that h=3.1 km.w.e and h=1.8 km.w.e, respectively. Figure [H 
shows muon energy spectra generated at Gran Sasso (black line) and at 
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Y2L (red line). Muon average energies are 276 GeV and 201 GeV for Gran 


Sasso and Y2L, which are consistent with the measured values [34 


3.2. Schematic layout of simulation geometry 

Inside the rock of a thickness of 20 m x 20 m x 20 m there is an air- 
filled cavern whose dimensions are 4mx4mx4m. A hemisphere with 
a radius of 10 m is assumed and generated muons are incident onto the 
shielded detector with an angle randomly selected according to an angular 
distribution proportional to cos 2 6 , as shown in Figure [2l The vertex of the 
muons is generated uniformly inside the square of 20 mx20 m area at the 
surface of the hemisphere. The size of the square is determined to cover 
the whole volume of the rock for all the angles of the generated muons. 
Figure [Tj shows the energy spectrum of generated muons at the hemisphere. 
The muon flux at Gran Sasso is considered as 2.7 x 10 -8 muons/cm 2 /s or 
2 x 10 7 muons /(20 m x 20 m)/2144 days and it is 2.7 x 10 - ' muons/cm 2 /s 
or 2 x 10' muons/(20 m x 20 m)/214 days at Y2L. 



Figure 2: Schematic layout of simulation geometry. 


3.3. Neutron production at the boundary of rock and cavern 

To validate the simulation for the neutron production we compared the 
neutron fluxes entering the cavern, which are produced in the rock. We also 
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included neutrons scattered back from the rock wall as well as shielding ma¬ 
terials surrounding the detector and reentering the cavern. Figure [3] shows 
the neutron energy spectra at the boundary between the rock and cavern, 
obtained from the simulation of muon propagation and interaction with ma¬ 
terials in the rock. There are four different shielding configurations inside the 
cavern, denoted as 1, 2, 3, and 4, defined in section I3~41 therefore, we can see 
the effect of the neutrons reflected back from the wall and shielding materials 
inside the cavern. Most are dominant below ~1 MeV and it depends on the 
outmost shield layer whether it is polyethylene (PE) or not. 

According to the simulation results neutron fluxes are 6.98xlO~ 10 n/cm 2 /s 
and 5.3 x 10” 9 n/cm 2 / s, for energies above 1 MeV, under shielding configu¬ 
ration 1 at Gran Sasso and under configuration 4 at Y2L, respectively. It is 
in good agreement with GEANT4 results of version 6.2 [35| and close to the 
FLUKA result [36] , 8.7 x 10 10 n/cm 2 /s, obtained from muon simulation in 
NaCl at Boulby. 



Neutron Energy (MeV) 


Figure 3: Neutron energy spectra at the boundary between the rock and cavern with dif¬ 
ferent shielding configurations inside the cavern. In the simulation, there are four shielding 
configurations: configuration 1 (thick black line), configuration 2 (thin black line), config¬ 
uration 3 (dotted blue line), and configuration 4 (thick red line). 
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3-4- Shielding configurations and shielding effects 

We assumed the Nal detectors used by DAMA/LIBRA, and therefore 
there are 5x5 segmented Nal crystals with dimensions of (10.2 x 10.2 x 
25.4) cm 3 for a single crystal. The target Nal crystals are inside a 10 cm 
copper shield and it is additionally surrounded by lead and polyethylene (PE) 
shields. 

To study the shielding effects, we set four different shielding configura¬ 
tions outside the copper shield: 

1. Lead (15 cm) + PE (30 cm) same as the DAMA setup 

2. PE (30 cm) + Lead (15 cm) 

3. Lead shield only (15 cm) 

4. PE (5 cm) + Lead (15 cm) + MD (30 cm) same as KIMS-Nal setup 

Shielding configuration 1 represents the DAMA setup at Gran Sasso 
depth, which uses 30 cm thick PE as the outmost shield layer. The order of 
lead and PE shield layers was switched for configuration 2. The shielding ef¬ 
fect was also tested with a lead shield only (configuration 3). Figure 0] shows 
shielded detector geometries for KIMS-Nal (left) and DAMA (right) setups. 
For the KIMS-Nal setup (configuration 4), 5 cm thick PE was put inside a 
lead shield and MD outside the lead shield, which is a muon detector filled 
with mineral oil (liquid paraffin). We expect that 30 cm thick MD can also 
shield neutrons as much as the PE shield did in the configuration 1. 



Figure 4: Shielded detector geometries for KIMS-Nal (left) and DAMA (right). 

However, neutrons can also be produced by muon interaction with shield¬ 
ing materials and should be considered in configuring shield layers. Fig¬ 
ure [5] (left) shows the primary vertex positions of neutrons in the x-y plane 
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of the shielding configuration 1, and it is found to be dense in the lead and 
copper shields. The right side of Figure 0 shows neutron energy spectra at 
boundaries between different shield layers. We defined three different bound¬ 
aries: rock/cavern (thick black line), outside copper shield (thin red line), 
and between copper and target crystals inside (dotted blue line). From the 
neutron spectra we found that the neutrons are built up in the lead and 
copper shields. 




Neutron Energy (MeV) 


Figure 5: Neutron vertex distribution (left) and neutron energy spectra at boundaries 
between different shield layers (right). 


To quantify the effects of neutrons with and without PE inside the lead 
shield we tested neutron energy spectra at boundaries with configuration 2, 
in which we changed the order of shield layers of PE and lead, as explained 
earlier. As shown in Figure [6] neutrons below ~10 MeV are reduced dramat¬ 
ically between PE and the copper shield. Therefore, it is more effective to 
shield neutrons below ~10 MeV by putting the PE shield inside the lead. 

In addition, we tested the shielding effects with respect to shielding con¬ 
figurations by comparing neutron energy spectra from configurations 1, 2, 
3, and 4 at two different boundaries. One boundary is outside the copper 
shield and another is inside the copper shield. Results at each boundary are 
shown in the left and the right of Figured There is little difference between 
configuration 1 and 3 for neutrons below ~0.1 MeV. If we consider the or¬ 
der of magnitude difference in muon fluxes at the depths of Gran Sasso and 
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Figure 6: Neutron energy spectra at boundaries between different shield layers with con¬ 
figuration 2. 



Neutron Energy (MeV) 


Neutron Energy (MeV) 


Figure 7: Effect of neutron shielding with different configurations. 


Y2L, the shielding effect with configuration 2 is very similar to the result 
with configuration 4, which used 5 cm thick PE inside the lead shield and 
30 cm thick MD instead of PE in configuration 1. As shown in the left of 
Figure [7] even 5 cm thick PE inside the lead is effective to shield neutrons 









































below ~1 MeV. Accordingly, 30 cm thick PE inside the lead shield (config¬ 
uration 2) is the most effective in shielding intermediate neutrons compared 
to other configurations. 

3.5. Neutron detection in Nal crystals 

From the simulation studies described in sections 13.31 and IXT1 we obtained 
the neutron flux at the boundaries of the shield layers and tested the shielding 
effect in each shielding material. We now consider the nuclear recoil event in 
Nal crystals due to neutron backgrounds. 



Neutron Energy Deposited in Nal crystal (MeV) 


Figure 8: Energy spectra deposited in Nal crystals. Results are obtained from the sim¬ 
ulation with shielding configuration f at Gran Sasso. The thick black line represents all 
events that deposited energy in the Nal crystal and the dotted black line represents single 
hit events. The green color indicates energy deposited by all nuclear recoil events, and the 
blue represents results from nuclear recoil single hit event only. 


Because we have 5x5 segmented Nal crystals an event can deposit its 
energy in two or more crystals. Therefore, we call an event a single hit event 
when it deposits the energy in a single Nal crystal and a multiple hit event 
when it exhibits energy deposition in two or more Nal crystals. Figure [8] 
shows the spectra of energy deposited in Nal crystals, which include not 
only nuclear recoil but also electromagnetic energy deposition. The thick 
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black line represents all events that deposited energy in the Nal crystal and 
the dotted black line represents single hit events. The green color represents 
energy deposited by all nuclear recoil events, and the blue color represents 
results for nuclear recoil single events only. The nuclear recoil events con¬ 
tribute only at the lowest energies. Figure [8] is obtained from the simulation 
with shielding configuration 1 at Gran Sasso, and for nuclear recoil single hit 
events it is much lower than the 10~ 4 cpd level for energies below 100 keV (in¬ 
set). 




Figure 9: Neutron energy deposited in Nal crystals by single hit events (left). The blue 
histogram is obtained from nuclear recoil single hit events, and the dotted black histogram 
results from single hit events, including electromagnetic energy deposition as well as nu¬ 
clear recoil. Quenching factors, 0.25 and 0.1, for Na and I are applied in the simulation. 
The correlation between neutron energy and nuclear recoil energy is shown as a scatter 
plot (right). 


The left of Figure [9] shows the spectra of visible energy deposited by all 
single hit events (dotted black line), using configuration 1. The thick blue 
line represents results from nuclear recoil single hit events only. A nuclear 
recoil event yields a lower number of scintillation photons than an electron 
recoil events for the same energy deposition. Measurements of the quenching 
factors for Na and I recoils in Nal(Tl) crystals have been reported and are ~25 
and ~10%, respectively [3sl]. Therefore, we applied 0.25 and 0.1 as quenching 
factors for Na and I in the simulation. The right graph in Figure [9] depicts 
the relationship between nuclear recoil energy and neutron kinetic energy as 
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a scatter plot, and the Na and I recoils are separated. 


4. Results 

Tabled] shows the summarized results for single hit event rates for the four 
shielding configurations. Single hit events include electromagnetic energy 
deposition as well as nuclear recoil energy deposition. The single hit event 
rates by the nuclear recoil event only are also obtained. 

The nuclear recoil single hit event rates at Y2L and Gran Sasso are 
(4.59±1.05) xl0~ 5 and (8.69±1.45) xlO -6 dru, respectively, and all single hit 
events rates are on the order of 10 -4 and 10 -5 dru for configurations 4 and 
1 . 


Tabic 1: Neutron background events in dru unit (events/day/kg/keV) in 5x5 segmented 
Nal crystals with different shielding configurations at two different underground sites 


Place 

Shielding 

configuration 

Single hit Events 
[/day/kg/keV] 

Nuclear Recoil 
Single hit Events 
[/day/kg/keV] 

Energy [keVee] 

2 - 10 

Energy [keVee] 

2 - 10 

GS 

1 

(1.62±0.20) xlO -5 

(8.69±1.45) xlO” 6 

GS 

2 

(6.52±1.25) xlO” 6 

(1.93±0.68) xlO -6 

GS 

3 

(1.84±0.21) xlO -5 

(1.35±0.18) xlO -5 

Y2L 

4 

(1.21±0.17) xlO -4 

(4.59±1.05) xlO -5 


5. Conclusion 

As shown in Table [T] it is very unlikely that the neutrons are the main 
source for the DAMA annual modulation, which was reported as (0.0112±0.0012) 
cpd/kg/keV in the (2-6) keV energy interval of the modulation amplitude 

m 

The DAMA group refuted Davis’s claim because the induced modulation 
amplitudes from neutrons induced by muons and by neutrinos are less than 
9 xlO” 6 cpd/kg/keV (2 xlO -5 cpd/kg/keV for neutrons produced in the 
lead shield) and less than 2 xlO -6 cpd/kg/keV, respectively from the sim¬ 
ulation, and they are less than 0.1% of the modulation amplitude measured 
by DAMA/LIBRA 0. 
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